function [pval_MHR_test, MAPE, maxSR, maxIR, PricingErrors] = MHR_test_linfac_for_dist(RF_sample, test_assets, M_factors, horizons)

% this code only uses excess return moments and an SDF M = 1 - b'(F-E[F])
% the p value is obtained using the main test given in proposition 2 in the
% paper, and as reported in Table 1. 

% the inputs are: 
% RF_sample (real net risk-free rate, i.e. a number like 0.005 for a 0.5% risk-free rate)
% test_assets (a T x N matrix of excess returns to N test assets, 5% should be 0.05)
% M_factors (a T x K matrix of excess returns to the factors in the SDF, 5% should be 0.05)
% horizons (horizons over which to jointly test)

% the outputs are:


K = length(M_factors(1,:));

% denote test assets as factors, potentially different from M_factors
factors = test_assets;

T = length(RF_sample);
Gross_Rf = 1 + RF_sample;
Gross_factors = Gross_Rf + factors;
Gross_all = cat(2,Gross_Rf, Gross_factors);

N = length(Gross_factors(1,:));
start_dt = max(horizons);
Ts = T - max(horizons) + 1;

% create inputs for SDF
EF = mean(M_factors(start_dt:end,:))';
VarF = (M_factors(start_dt:end,:)-EF')'*(M_factors(start_dt:end,:)-EF')/Ts;

% constant param M
b0 = VarF \ EF;

% SDF that prices single-horizon returns to the M_factors 
M = 1 - (M_factors - EF')* b0;

% Construct moments
H = length(horizons);

% Construct MHR moments. Note that we are not testing for pricing of the
% risk-free asset here as the linear factor models we test in the paper are
% designed to price excess returns.
mH = horizons(end)-1;
zt = zeros(Ts,N,mH);
for jh = 1:mH
    for jn = 1:N
        zt(:,jn,jh) = GeoSum( M(start_dt-jh:end-1).*Gross_all(start_dt-jh:end-1,jn+1),jh);
    end
end

% now, create covariance conditions, as well as initial conditions
ft_new = zeros(Ts,N,H);
for jh = 1:H
    for jn = 1:N
        if horizons(jh) == 1
            ft_new(:,jn,jh) = 1 .* ( M(start_dt:end).*factors(start_dt:end,jn) );
        else
            ft_new(:,jn,jh) = sum(squeeze(zt(:,jn,horizons(1):horizons(jh)-1)),2) .*...
                ( M(start_dt:end).*factors(start_dt:end,jn) );
        end
    end
end

ft = reshape(ft_new,Ts,N*H);

% add b'mu moment
ft = cat(2,(M_factors(start_dt:end,:) - EF') ,ft);

gT = mean(ft)';

% K factor means and K b0 parameters estimated means -2*K dof
deg_free = length(gT) - 2*K;

% re-size gT for later alpha code, remove the first K moments that estimate
% factor means
gT = gT(K+1:end);

UMVE = M_factors(start_dt:end,:) * b0;
maxSR = sqrt(12) * mean(UMVE) / std(UMVE);

% get implied alpha and sigma inverse simply by running GRS type
% regressions as these are equivalent in SHR setting.

% get alphas:
alphas = gT(N+1:end) / mean(M(start_dt:end));

% use GRS/BJS regressions to get alphas and residual covariance matrix
% Test assets
Test_ExRet = zeros(Ts,N,H-1);
Inst_zt = zeros(Ts,N,H-1);
alpha = zeros(N,H-1);
eps = zeros(Ts,N,H-1);
M_times_eps = zeros(Ts,N,H-1);
R2s = zeros(N,H-1);
for jn = 2:N+1
    for jh = 2:H
        Inst_zt(:,jn-1,jh-1) = sum(squeeze(zt(:,jn-1,horizons(1):horizons(jh)-1)),2);
        Test_ExRet(:,jn-1,jh-1) = Inst_zt(:,jn-1,jh-1) .* (Gross_factors(start_dt:end,jn-1) -  Gross_Rf(start_dt:end,1));
        X = cat(2,ones(Ts,1),M_factors(start_dt:end,:));
        Y = Test_ExRet(:,jn-1,jh-1);
        beta_temp = (X'*X) \ X'*Y;
        eps_temp = Y - X * beta_temp;
        alpha(jn-1,jh-1) = beta_temp(1);
        R2s(jn-1,jh-1) = var(eps_temp)/var(Y);
        eps(:,jn-1,jh-1) = eps_temp;
        
        M_times_eps(:,jn-1,jh-1) = eps_temp .* M(start_dt:end,1);

    end
end
temp_ret = reshape(Test_ExRet, Ts, (H-1)*N,1);

% data for figure(1)
PricingErrors = nan(H-1,N);
for jn = 1:N
    PricingErrors(:,jn) = 12*mean(M(start_dt:end))*...
        alpha(jn,:) ./ horizons(2:end,1)';
end

M_times_eps = reshape(M_times_eps,Ts,N*(H-1),1);
eps = reshape(eps,Ts,N*(H-1),1);
alpha = reshape(alpha,N*(H-1),1);
V_M_times_eps = M_times_eps'*M_times_eps / Ts;
chi2_stat = Ts * (alpha' / V_M_times_eps) * alpha;
pval_MHR_test = 1-chi2cdf(chi2_stat,deg_free);


% test that alphas from GMM are the same as those from GRS regressions
if sum(abs(alphas-alpha)) > 1e-8
    disp('!!!!!!!!!!!!!Sum of GMM vs GRS alphas too large!!!!!!!!!');
end

% max IR^2
eps = reshape(eps,Ts,N*(H-1),1);
Sigma = cov(eps);
maxIR = sqrt( alpha' * inv(Sigma) * alpha * 12);

% mean pricing errors
MAPE = mean(mean(abs(PricingErrors)));


